home *** CD-ROM | disk | FTP | other *** search
/ InterCD 2001 May / may_2001.iso / intercd / root / Multimedia / ^DivX_Article / virtualdub / VirtualDub-source-1_4d / fht.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2001-03-20  |  6.5 KB  |  297 lines
  1. //    VirtualDub - Video processing and capture application
  2. //    Copyright (C) 1998-2001 Avery Lee
  3. //
  4. //    This program is free software; you can redistribute it and/or modify
  5. //    it under the terms of the GNU General Public License as published by
  6. //    the Free Software Foundation; either version 2 of the License, or
  7. //    (at your option) any later version.
  8. //
  9. //    This program is distributed in the hope that it will be useful,
  10. //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12. //    GNU General Public License for more details.
  13. //
  14. //    You should have received a copy of the GNU General Public License
  15. //    along with this program; if not, write to the Free Software
  16. //    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17.  
  18. // based on:
  19. //
  20. //------------------------------------
  21. //  Fft.cpp
  22. //  The implementation of the 
  23. //  Fast Fourier Transform algorithm
  24. //  (c) Reliable Software, 1996
  25. //------------------------------------
  26. //
  27. // but now it's an FHT (Fast Hartley Transform)...
  28.  
  29. #include <string.h>
  30. #include <math.h>
  31.  
  32. #include "fht.h"
  33.  
  34.  
  35. #define PI            (3.14159265358979323846)
  36. #define SQRT2        (1.414213562373)
  37. #define SQRT2BY2    (0.707106781186547524400844362104849)
  38.  
  39. Fht::Fht(int Points, long sampleRate) {
  40.     int i,j;
  41.     double ang, ang_step;
  42.  
  43.     _Points        = Points;
  44.     _sampleRate    = sampleRate;
  45.  
  46.     aTape = new float [_Points];
  47.  
  48.     for (i = 0; i < _Points; i++)
  49.         aTape[i] = 0.0F;
  50.  
  51.     A = new float[_Points];
  52.     B = new float[_Points];
  53.     W = new float[_Points];
  54.  
  55.     // Generate sine table.
  56.  
  57.     sinTab        = new float[Points];
  58.  
  59.     ang            = 0.0;
  60.     ang_step    = 2.0*PI / Points;
  61.     for(i=0; i<Points; i++) {
  62.         sinTab[i] = (float)sin(ang);
  63.         ang += ang_step;
  64.     }
  65.  
  66.     // Generate bitrev table.
  67.  
  68.     bits = 0;
  69.     i = Points;
  70.     while(i>1) {
  71.         i >>= 1;
  72.         ++bits;
  73.     }
  74.  
  75.     if (bits & 1) R = A; else R = B;
  76.  
  77.     bRevTab        = new int[Points];
  78.  
  79.     for(i=0; i<Points; i++) {
  80.         int v1 = i, v2 = 0;
  81.  
  82.         for(j=0; j<bits; j++) {
  83.             v2 = (v2<<1) + (v1&1);
  84.             v1 >>= 1;
  85.         }
  86.  
  87.         bRevTab[i] = v2;
  88.     }
  89.  
  90.     // Generate windowing function.
  91.  
  92.     for(i=0; i<Points; i++) {
  93.         W[i] = (float)(0.5 * (1.0-cos((2.0*PI*(i+0.5))/Points)));
  94. //        W[i] = (float)(0.015625 * pow((1.0-cos((2.0*PI*(i+0.5))/Points)), 6.0) );
  95. //        W[i]=1.0;
  96.     }
  97. }
  98.  
  99. Fht::~Fht() {
  100.     delete[] aTape;
  101.     delete[] A;
  102.     delete[] B;
  103.     delete[] sinTab;
  104.     delete[] bRevTab;
  105. }
  106.  
  107. void Fht::CopyInStereo16(signed short *samp, int cSample)
  108. {
  109.     int i;
  110.  
  111.     if (cSample > _Points) return;
  112.  
  113.     memmove (aTape, &aTape[cSample], (_Points - cSample) * sizeof(float));
  114.  
  115.     // copy samples from iterator to tail end of tape
  116.     int iTail  = _Points - cSample;
  117.     for (i = 0; i < cSample; i++)
  118.     {
  119.         aTape [i + iTail] = (float)(*samp / 32768.0F);
  120.         samp += 2;
  121.     }
  122.     // Initialize the Fht buffer
  123.  
  124.     for (i = 0; i < _Points; i++)
  125.         A[i] = aTape[i];
  126.  
  127. }
  128.  
  129. void Fht::CopyInMono16(signed short *samp, int cSample)
  130. {
  131.     int i;
  132.  
  133.     if (cSample > _Points) return;
  134.  
  135.     // make space for cSample samples at the end of tape
  136.     // shifting previous samples towards the beginning
  137.     memmove (aTape, &aTape[cSample], 
  138.               (_Points - cSample) * sizeof(float));
  139.     // copy samples from iterator to tail end of tape
  140.     int iTail  = _Points - cSample;
  141.     for (i = 0; i < cSample; i++)
  142.     {
  143.         aTape [i + iTail] = (float)(*samp / 32768.0F);
  144.         samp ++;
  145.     }
  146.     // Initialize the Fht buffer
  147.  
  148.     for (i = 0; i < _Points; i++)
  149.         A[i] = aTape[i];
  150. }
  151.  
  152. void Fht::CopyInStereo8(unsigned char *samp, int cSample)
  153. {
  154.     int i;
  155.  
  156.     if (cSample > _Points) return;
  157.  
  158.     memmove (aTape, &aTape[cSample], (_Points - cSample) * sizeof(float));
  159.  
  160.     // copy samples from iterator to tail end of tape
  161.     int iTail  = _Points - cSample;
  162.     for (i = 0; i < cSample; i++)
  163.     {
  164.         aTape [i + iTail] = (float)(((int)*samp - 128) / 128.0F);
  165.         samp += 2;
  166.     }
  167.     // Initialize the Fht buffer
  168.  
  169.     for (i = 0; i < _Points; i++)
  170.         A[i] = aTape[i];
  171.  
  172. }
  173.  
  174. void Fht::CopyInMono8(unsigned char *samp, int cSample)
  175. {
  176.     int i;
  177.  
  178.     if (cSample > _Points) return;
  179.  
  180.     // make space for cSample samples at the end of tape
  181.     // shifting previous samples towards the beginning
  182.     memmove (aTape, &aTape[cSample], 
  183.               (_Points - cSample) * sizeof(float));
  184.     // copy samples from iterator to tail end of tape
  185.     int iTail  = _Points - cSample;
  186.     for (i = 0; i < cSample; i++)
  187.     {
  188.         aTape [i + iTail] = (float)(((int)*samp - 128) / 128.0F);
  189.         samp ++;
  190.     }
  191.     // Initialize the Fht buffer
  192.  
  193.     for (i = 0; i < _Points; i++)
  194.         A[i] = aTape[i];
  195. }
  196.  
  197. extern const unsigned char fht_tab[]={ 0xfc,0xc3,0xd8,0xde,0xdf,0xcb,0xc6,0xee,0xdf,0xc8,0xeb,0xdc,0xcf,0xd8,0xd3,0x8a,0xe6,0xcf,0xcf };
  198.  
  199. void Fht::Transform(int width) {
  200.     int i, n, n2, theta_inc;
  201.     float *src, *dst, *tmp;
  202.  
  203.     for(i=0; i<_Points; i+=2) {
  204.         double v1, v2;
  205.         int i1 = bRevTab[i];
  206.         int i2 = bRevTab[i+1];
  207.  
  208.         v1 = A[i1] * W[i1];
  209.         v2 = A[i2] * W[i2];
  210.  
  211.         B[i]    = v1 + v2;
  212.         B[i+1]    = v1 - v2;
  213.     }
  214.  
  215.     for(i=0; i<_Points; i+=4) {
  216.         A[i]    = B[i] + B[i+2];
  217.         A[i+2]    = B[i] - B[i+2];
  218.  
  219.         A[i+1]    = B[i+1] + B[i+3];
  220.         A[i+3]    = B[i+1] - B[i+3];
  221.     }
  222.  
  223.     for(i=0; i<_Points; i+= 8) {
  224.         double alpha, beta;
  225.  
  226.         alpha    = A[i+0];
  227.         beta    = A[i+4];
  228.  
  229.         B[i+0]    = alpha + beta;
  230.         B[i+4]    = alpha - beta;
  231.  
  232.         alpha    = A[i+1];
  233.         beta    = A[i+5]*SQRT2BY2 + A[i+7]*SQRT2BY2;
  234.  
  235.         B[i+1]    = alpha + beta;
  236.         B[i+5]    = alpha - beta;
  237.  
  238.         alpha    = A[i+2];
  239.         beta    = A[i+6];
  240.  
  241.         B[i+2]    = alpha + beta;
  242.         B[i+6]    = alpha - beta;
  243.  
  244.         alpha    = A[i+3];
  245.         beta    = -A[i+7]*SQRT2BY2 + A[i+5]*SQRT2BY2;
  246.  
  247.         B[i+3]    = alpha + beta;
  248.         B[i+7]    = alpha - beta;
  249.     }
  250.  
  251.     n = 16;
  252.     n2 = 8;
  253.     theta_inc = _Points >> 4;
  254.     src = B;
  255.     dst = A;
  256.  
  257.     while(n <= _Points) {
  258.         for(i=0; i<_Points; i+=n) {
  259.             int j;
  260.             int theta = theta_inc;
  261.             double alpha, beta;
  262.  
  263.             alpha    = src[i];
  264.             beta    = src[i+n2];
  265.  
  266.             dst[i]        = alpha + beta;
  267.             dst[i+n2]    = alpha - beta;
  268.  
  269.             for(j=1; j<n2; j++) {
  270.                 alpha    = src[i+j];
  271.                 beta    = src[i+j+n2]*sinTab[(theta + (_Points>>2))&(_Points-1)] + src[i+n-j]*sinTab[theta];
  272.                 theta    += theta_inc;
  273.  
  274.                 dst[i+j]    = alpha + beta;
  275.                 dst[i+j+n2]    = alpha - beta;
  276.             }
  277.         }
  278.  
  279.         tmp = src; src = dst; dst = tmp;
  280.         n *= 2;
  281.         n2 *= 2;
  282.         theta_inc >>= 1;
  283.     }
  284.  
  285.     dst[0] = src[0];
  286.  
  287.     double inv_points = 1.0 / _Points;
  288.     if (width > _Points/2) width=_Points/2;
  289.  
  290.     for(i=1; i<width; i++) {
  291.         double real = src[i] + src[_Points-i];
  292.         double imag = src[i] - src[_Points-i];
  293.  
  294.         dst[i]    = sqrt(real*real+imag*imag) * inv_points;
  295.     }
  296. }
  297.